1 Objective

The objective of this notebook is to perform pre-processing and dimensional reduction on both assays independently, using standard approaches for RNA and ATAC-seq data. Then, we will follow the “Joint RNA and ATAC analysis: 10x multiomic” vignette from Signac to obtain a joint visualization using both modalities.

2 Pre-processing

2.1 Load packages

library(Signac)
library(Seurat)
library(GenomicRanges)
library(future)
library(SeuratWrappers)
library(harmony)
library(EnsDb.Hsapiens.v86)
library(stringr)
library(dplyr)
library(ggplot2)

set.seed(173)

2.2 Parameters

path_to_obj <- here::here("multiome/results/R_objects/4.tonsil_multiome_filtered_combined_with_metadata.rds")
path_to_save <- here::here("multiome/results/R_objects/5.tonsil_multiome_integrated_using_wnn.rds")

2.3 Load Multiome filtered data

tonsil_filtered <- readRDS(path_to_obj)
tonsil_filtered
## An object of class Seurat 
## 188447 features across 47150 samples within 2 assays 
## Active assay: peaks (151846 features, 0 variable features)
##  1 other assay present: RNA

3 Without Harmony integration

3.1 scATAC

3.1.1 Normalization and linear dimensional reduction

DefaultAssay(tonsil_filtered) <- "peaks"
tonsil_filtered <- RunTFIDF(tonsil_filtered)
tonsil_filtered <- FindTopFeatures(tonsil_filtered, min.cutoff = "q0")
tonsil_filtered <- RunSVD(tonsil_filtered)
DepthCor(tonsil_filtered)

3.1.2 UMAP representation

tonsil_filtered <- RunUMAP(
  tonsil_filtered,
  dims = 2:40,
  reduction = "lsi",
  reduction.name = "umap.atac",
  reduction.key = "atacUMAP_"
)
DimPlot(
  tonsil_filtered,
  reduction = "umap.atac",
  group.by = "library_name",
  pt.size = 0.1
) + NoLegend()

3.2 scRNA

3.2.1 Normalization and linear dimensional reduction

DefaultAssay(tonsil_filtered) <- "RNA"
tonsil_filtered <- NormalizeData(
  tonsil_filtered,
  normalization.method = "LogNormalize",
  scale.factor = 1e4
)
tonsil_filtered <- tonsil_filtered %>%
  FindVariableFeatures(nfeatures = 3000) %>%
  ScaleData() %>% 
  RunPCA() 

3.2.2 UMAP representation

tonsil_filtered <- RunUMAP(
  tonsil_filtered,
  dims = 1:30,
  reduction = "pca",
  reduction.name = "umap.rna",
  reduction.key = "rnaUMAP_"
)
DimPlot(
  tonsil_filtered,
  reduction = "umap.rna",
  group.by = "library_name",
  pt.size = 0.1
) + NoLegend()

3.3 Joint

 #tonsil_filtered <- FindMultiModalNeighbors(
  #  tonsil_filtered,
  #  reduction.list = list("pca", "lsi"),
  #  dims.list = list(1:30, 2:40)
 #)
#  tonsil_filtered <- RunUMAP(
#    tonsil_filtered,
#    nn.name = "weighted.nn",
#    reduction.name = "wnn.umap",
#    reduction.key = "wnnUMAP_"
#  )
#  DimPlot(
 #   tonsil_filtered,
#    reduction = "wnn.umap",
#    group.by = "library_name",
#    pt.size = 0.1
#  ) + NoLegend()

4 With Harmony integration

4.1 scATAC

DefaultAssay(tonsil_filtered) <- "peaks"
tonsil_filtered <- RunHarmony(
  object = tonsil_filtered,
  reduction = "lsi",
  dims = 2:40,
  group.by.vars = "gem_id",
  assay.use = "peaks",
  project.dim = FALSE,
  reduction.save = "harmony_peaks"
)

4.1.1 UMAP representation

tonsil_filtered <- RunUMAP(
  tonsil_filtered,
  dims = 2:40,
  reduction = "harmony_peaks",
  reduction.name = "umap.atac",
  reduction.key = "atacUMAP_"
)
DimPlot(
  tonsil_filtered,
  reduction = "umap.atac",
  group.by = "library_name",
  pt.size = 0.1
) + NoLegend()

4.2 scRNA

DefaultAssay(tonsil_filtered) <- "RNA"
tonsil_filtered <- RunHarmony(
  object = tonsil_filtered,
  reduction = "pca",
  dims = 1:30,
  group.by.vars = "gem_id",
  assay.use = "RNA",
  project.dim = FALSE,
  reduction.save = "harmony_RNA"
)
tonsil_filtered <- RunUMAP(
  tonsil_filtered,
  dims = 1:30,
  reduction = "harmony_RNA",
  reduction.name = "umap.rna",
  reduction.key = "rnaUMAP_"
)
DimPlot(
  tonsil_filtered,
  reduction = "umap.rna",
  group.by = "library_name",
  pt.size = 0.1
) + NoLegend()

4.3 Joint

#  tonsil_filtered <- FindMultiModalNeighbors(
#    tonsil_filtered,
 #   reduction.list = list("harmony_peaks", "harmony_RNA"),
#    dims.list = list(2:40, 1:30)
#  )
#  tonsil_filtered <- RunUMAP(
#    tonsil_filtered,
#    nn.name = "weighted.nn",
#    reduction.name = "wnn.umap",
#    reduction.key = "wnnUMAP_"
#  )
#  DimPlot(
#    tonsil_filtered,
#    reduction = "wnn.umap",
#    group.by = "library_name",
#    pt.size = 0.1
#  ) + NoLegend()

5 Save

We will save the resulting object and use it in the following notebook to exclude doublets:

saveRDS(tonsil_filtered, path_to_save)

6 Session Information

sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: macOS Big Sur 10.16
## 
## Matrix products: default
## BLAS/LAPACK: /Users/pauli/opt/anaconda3/envs/Tonsil_atlas/lib/libopenblasp-r0.3.10.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] ggplot2_3.3.2             dplyr_1.0.2               stringr_1.4.0             EnsDb.Hsapiens.v86_2.99.0 ensembldb_2.12.1          AnnotationFilter_1.12.0   GenomicFeatures_1.40.1    AnnotationDbi_1.50.3      Biobase_2.48.0            harmony_1.0               Rcpp_1.0.5                SeuratWrappers_0.3.0      future_1.20.1             GenomicRanges_1.40.0      GenomeInfoDb_1.24.0       IRanges_2.22.1            S4Vectors_0.26.0          BiocGenerics_0.34.0       Seurat_3.9.9.9010         Signac_1.1.0.9000         BiocStyle_2.16.1         
## 
## loaded via a namespace (and not attached):
##   [1] reticulate_1.18             tidyselect_1.1.0            RSQLite_2.2.1               htmlwidgets_1.5.2           grid_4.0.3                  BiocParallel_1.22.0         Rtsne_0.15                  munsell_0.5.0               codetools_0.2-17            ica_1.0-2                   miniUI_0.1.1.1              withr_2.3.0                 colorspace_2.0-0            OrganismDbi_1.30.0          knitr_1.30                  rstudioapi_0.12             ROCR_1.0-11                 tensor_1.5                  listenv_0.8.0               labeling_0.4.2              GenomeInfoDbData_1.2.3      polyclip_1.10-0             farver_2.0.3                bit64_4.0.5                 rprojroot_2.0.2             parallelly_1.21.0           vctrs_0.3.4                 generics_0.1.0              xfun_0.18                   biovizBase_1.36.0           BiocFileCache_1.12.1        lsa_0.73.2                  ggseqlogo_0.1               R6_2.5.0                    rsvd_1.0.3                  bitops_1.0-6                spatstat.utils_1.17-0       reshape_0.8.8               DelayedArray_0.14.0         assertthat_0.2.1            promises_1.1.1              scales_1.1.1               
##  [43] nnet_7.3-14                 gtable_0.3.0                globals_0.13.1              goftest_1.2-2               ggbio_1.36.0                rlang_0.4.8                 RcppRoll_0.3.0              splines_4.0.3               rtracklayer_1.48.0          lazyeval_0.2.2              dichromat_2.0-0             checkmate_2.0.0             BiocManager_1.30.10         yaml_2.2.1                  reshape2_1.4.4              abind_1.4-5                 backports_1.2.0             httpuv_1.5.4                Hmisc_4.4-1                 RBGL_1.64.0                 tools_4.0.3                 bookdown_0.21               ellipsis_0.3.1              RColorBrewer_1.1-2          ggridges_0.5.2              plyr_1.8.6                  base64enc_0.1-3             progress_1.2.2              zlibbioc_1.34.0             purrr_0.3.4                 RCurl_1.98-1.2              prettyunits_1.1.1           rpart_4.1-15                openssl_1.4.3               deldir_0.2-3                pbapply_1.4-3               cowplot_1.1.0               zoo_1.8-8                   SummarizedExperiment_1.18.1 ggrepel_0.8.2               cluster_2.1.0               here_1.0.1                 
##  [85] magrittr_1.5                RSpectra_0.16-0             data.table_1.13.2           lmtest_0.9-38               RANN_2.6.1                  SnowballC_0.7.0             ProtGenerics_1.20.0         fitdistrplus_1.1-1          matrixStats_0.57.0          hms_0.5.3                   patchwork_1.1.0             mime_0.9                    evaluate_0.14               xtable_1.8-4                XML_3.99-0.3                jpeg_0.1-8.1                gridExtra_2.3               compiler_4.0.3              biomaRt_2.44.4              tibble_3.0.4                KernSmooth_2.23-17          crayon_1.3.4                htmltools_0.5.0             mgcv_1.8-33                 later_1.1.0.1               Formula_1.2-4               tidyr_1.1.2                 DBI_1.1.0                   tweenr_1.0.1                dbplyr_1.4.4                MASS_7.3-53                 rappdirs_0.3.1              Matrix_1.2-18               igraph_1.2.6                pkgconfig_2.0.3             GenomicAlignments_1.24.0    foreign_0.8-80              plotly_4.9.2.1              xml2_1.3.2                  XVector_0.28.0              VariantAnnotation_1.34.0    digest_0.6.27              
## [127] sctransform_0.3.1           RcppAnnoy_0.0.16            graph_1.66.0                spatstat.data_1.4-3         Biostrings_2.56.0           rmarkdown_2.5               leiden_0.3.5                fastmatch_1.1-0             htmlTable_2.1.0             uwot_0.1.8.9001             curl_4.3                    shiny_1.5.0                 Rsamtools_2.4.0             lifecycle_0.2.0             nlme_3.1-150                jsonlite_1.7.1              viridisLite_0.3.0           askpass_1.1                 BSgenome_1.56.0             pillar_1.4.6                lattice_0.20-41             GGally_2.0.0                fastmap_1.0.1               httr_1.4.2                  survival_3.2-7              remotes_2.2.0               glue_1.4.2                  spatstat_1.64-1             png_0.1-7                   bit_4.0.4                   ggforce_0.3.2               stringi_1.5.3               blob_1.2.1                  latticeExtra_0.6-29         memoise_1.1.0               irlba_2.3.3                 future.apply_1.6.0